Librerías

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(readxl)
library(plotly)
## Warning: package 'plotly' was built under R version 3.5.1
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.1
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

CAPITULO 1

Mezclas

Al tener una coleccion de datos observable, y estudiarlo a traves de correlaciones, podemos encontrar que puede tener 1 o mas parametros que rigen el comportamiento de los datos. La herramienta que nos permite el analisis de los parametros y como estos impactan el comportamiento de los datos es Hiden Markiv models o H.M.M.

Los H.M.M son modelos de distribucion probabilisticos de proposito general para una o más variables. Y es especifico para casos de sieres de tiempo, categoricas y continuas.

se debe evaluar cual es la distribucion que mejor se apega a los datos, y para el proposito de este documento (estudio de terremotos) se utilizara la distribución de poisson.

A continuacion un de la cuenta de sismos por anio, estos datos nos permitiran ejemplificar las mezclas y como la necesidad de los datos para metricos nos insta a utilizar H.M.M.

Año Sismos por año
1900 13
1901 14
1902 10
1903 9
1904 13
1905 18
1906 24
1907 18
1908 11
1909 19
1910 28
1911 18
1912 14
1913 18
1914 17
1915 17
1916 25
1917 18
1918 17
1919 13
1920 8
1921 12
1922 14
1923 22
1924 17
1925 15
1926 16
1927 17
1928 16
1929 16
1930 8
1931 23
1932 12
1933 13
1934 20
1935 20
1936 16
1937 20
1938 22
1939 17
1940 18
1941 20
1942 17
1943 32
1944 23
1945 14
1946 23
1947 14
1948 17
1949 16
1950 24
1951 13
1952 12
1953 15
1954 9
1955 12
1956 10
1957 23
1958 9
1959 12
1960 14
1961 11
1962 9
1963 17
1964 12
1965 13
1966 11
1967 10
1968 23
1969 15
1970 21
1971 18
1972 15
1973 10
1974 16
1975 14
1976 17
1977 14
1978 17
1979 14
1980 13
1981 11
1982 10
1983 15
1984 13
1985 13
1986 11
1987 14
1988 10
1989 6
1990 19
1991 13
1992 13
1993 12
1994 17
1995 20
1996 14
1997 17
1998 12
1999 18
2000 12
2001 14
2002 13
2003 14
2004 15
2005 11
2006 11
2007 14 
library(dplyr)
library(plotly)
datos <- data.frame(anio =c(1900,1901,1902,1903,1904,1905,1906,1907,1908,1909,1910,1911,1912,1913,1914,1915,1916,1917,1918,1919,1920,1921,1922,1923,1924,1925,1926,1927,1928,1929,1930,1931,1932,1933,1934,1935,1936,1937,1938,1939,1940,1941,1942,1943,1944,1945,1946,1947,1948,1949,1950,1951,1952,1953,1954,1955,1956,1957,1958,1959,1960,1961,1962,1963,1964,1965,1966,1967,1968,1969,1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,
1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007), frecuencia = c(13,14,10,9,13,18,24,18,11,19,28,18,14,18,17,17,25,18,17,13,8,12,14,22,17,15,16,17,16,16,8,23,12,13,20,20,16,20,22,17,18,20,17,32,23,14,23,14,17,16,24,13,12,15,9,12,10,23,9,12,14,11,9,17,12,13,11,10,23,15,21,18,15,10,16,14,17,14,17,14,13,11,10,15,13,13,11,14,10,6,19,13,13,12,17,20,14,17,12,18,12,14,13,14,15,11,11,14))

plot_ly(data = datos, x=~anio,  y=~frecuencia, type = 'scatter', mode = 'lines+markers')
plot_ly(data = datos, x=~anio,  y=~frecuencia, type = 'bar')

Modelo de mezclas independientes

Al utilizar la distribucion de poisson \(p(x) = \frac{exp(-\lambda)\lambda^x}{x!}\) tenemos que la media es de 15, pero segun poisson la media y la varianza son iguales, para este ejemplo la varianza es de 20.33, no coinciden con la distribucion de Pisson, sin embargo el fenomeno satisface las condiciones para utilizar la distribucion de Poisson. Se plantea la hipostesis, que el resultado esta siendo afectado por 1 de 2 distribuciones de Poisson con medias \(\lambda_1, \lambda_2\).

La seleccion de que media impacta el proceso esta dada por un mecanismo aleatorio, suponemos que la media \(\lambda_1\) es seleccionada con probabilidad \(\delta_1\) y \(\lambda_2\) con probabilidad \(\delta_2= 1-\delta_1\).

Vamos a definir un mecanismo C para definir las mezclas donde:

\(C=\left\{ \begin{array}{lcc} 1 & con & probabilidad &\delta_1 \\ \\ 2 & con & probabilidad &\delta_2 \\ \end{array} \right.\)

La funcion de probabilidad para \(X\) esta dada por \[\sum_{i=1}{\delta_i p_i}\]

Estimacion del parametro

Para la estimacion de los parametros un de una distribucion mixta esta dada por la maximo parentesco (Maximum Likeleyhood) la cual es la misma es para casos discretos o continuos. \[L(\theta_1,...,\theta_m,\delta_1,...,\delta_m)= \prod_{j=1}^{n}\sum_{i=1}^{m}\delta_i p_i(x_j,\theta_i)\]. Los \(\theta\) son los vectores de los parametros de la distribucion y los \(\delta\) son los parametros de mezcla, los cuales suman 1; los \(x_i,...,x_n\) son las observaciones.

Para el caso de m=2 (2 componentes de la distribucion de Poisson) con medias \(\lambda_1\) y \(\lambda_2\), y sean \(\delta_1\) y \(\delta_2=1-\delta_1\). Tenemos entonces que la distribucion p es: \[p(x)=\delta_1 \frac{\lambda_i^{x_i} exp(-\lambda_i)}{x!}+(1-\delta_1)\frac{\lambda_2^{x_i} exp(-\lambda_2)}{x_i!}\]

Asumiendo que \(\delta_2 = 1-\delta_1\), los parametros a buscar son \(\lambda_1, \lambda_2,\delta_1\) y el maximo parentesco es:\[L(\lambda_1, \lambda_2, \delta_1 | x_1,...,x_n) = \prod_{i=1}^n (\delta_1 \frac{\lambda_1^{x_i}e^{-\lambda_1}}{x_i!} + (1-\delta_1) \frac{\lambda_2^{x_i}e^{-\lambda_2}}{x_i!})\]

mllk <- function(wpar,x){
  zzz <- w2n(wpar)
  -sum(log(outer(x,zzz$lambda,dpois)%*%zzz$delta))
}

n2w <- function(lambda,delta){
  log(c(lambda, delta[-1]/(1-sum(delta[-1]))))
}

w2n <- function(wpar){
  m<- (length(wpar)+1)/2
  lambda<- exp(wpar[1:m])
  delta <- exp(c(0,wpar[(m+1):(2*m-1)]))
  return (list(lambda=lambda,delta=delta/sum(delta)))
}

likelihood <- function(parms){
  return(prod(outer(x,params$lambda,dpois) %*% params$delta))
}

Veamos con \(\delta_1\) y \(\delta_2\)

x <- datos$frecuencia
wpar <- n2w(c(10,20),c(0.5,0.5))
params<-w2n(nlm(mllk,wpar,x)$estimate)
## Warning in nlm(mllk, wpar, x): NA/Inf replaced by maximum positive value
L = likelihood(params)
params <-  c(params, likelihood=L, "log L"=-log(L))
as.data.frame(params)
##     lambda      delta    likelihood    log.L
## 1 22.59287 0.09397933 2.096784e-135 310.1086
## 2 14.59052 0.90602067 2.096784e-135 310.1086

Veamos con \(\delta_1\), \(\delta_2\) y \(\delta_3\)

x <- datos$frecuencia
wpar <- n2w(c(10,20,30),c(1,1,1)/3)
params<-w2n(nlm(mllk,wpar,x)$estimate)
## Warning in nlm(mllk, wpar, x): NA/Inf replaced by maximum positive value
L = likelihood(params)
params <-  c(params, likelihood=L, "log L"=-log(L))
as.data.frame(params)
##     lambda      delta    likelihood    log.L
## 1 14.59060 0.31459256 2.096784e-135 310.1086
## 2 14.59050 0.59143019 2.096784e-135 310.1086
## 3 22.59295 0.09397725 2.096784e-135 310.1086

Veamos con \(\delta_1\), \(\delta_2\), \(\delta_3\) y \(\delta_4\)

x <- datos$frecuencia
wpar <- n2w(c(10,20,35,40),c(1,1,1,1)/4)
params<-w2n(nlm(mllk,wpar,x)$estimate)
L = likelihood(params)
params <-  c(params, likelihood=L, "log L"=-log(L))
as.data.frame(params)
##     lambda      delta    likelihood    log.L
## 1 14.59058 0.25197174 2.096784e-135 310.1086
## 2 14.59051 0.47784478 2.096784e-135 310.1086
## 3 14.59061 0.17620764 2.096784e-135 310.1086
## 4 22.59299 0.09397583 2.096784e-135 310.1086

Como podemos apreciar a partir del \(\delta_2\) no hay ganancia, y el algoritmo esta descomponiendo la \(\delta_1\) para llenar las variables faltantes, por ello -log(L) da el mismo resultado.

Una mejor apreciacion de mejor ajuste es a traves de la esperanza matematica y la varianza \(E(X)=\sum_i\delta_i\sigma_i\) \(Var(X)=E(X^2)-(E(X))^2\) siendo \(E(X^2)=\sum_i\delta_i(\lambda_i+\lambda_i^2)\)

Cadenas de Markov

Una secuencia de variales aleatorias discretas \({C_t:t\in \mathbb{N}}\) son consideradas como una cadena de Markov si para todos los \(t\in \mathbb{N}\) se satisface \[Pr(C_{t+1}|C_t...C_1)=Pr(C_{t+1}|C_t)\]

Es decir, el estado actual depende solo y unicamente del estado anterior, o en otras palabras; no tiene memoria. Por razones de espacio podemos escribir \(C^{(t)}\) como la historio \((C_1,C_2,...,C_t)\) o sea \(Pr(C_{t+1}|C^{t})=Pr(C_{t+1}|C_t)\)

Las cantidades asociadas a los estados de las cadenas de Markov son las probailidades condicionadas tambien llamdas probabilidades de transicion. \(Pr(C_{s+t}=j|C_s=i)\)

Si las probabilidades no dependen de \(s\) entonces la matriz es llamada “homogenea”, de lo contrario es “no homogenea”. Asumiremos que las matrices son homogeneas a menos que se indique lo contrario. \(\gamma_{ij}=Pr(C_{s+t}=j|C_s=i)\)

La matriz \(\Gamma(t)\) es definida como la matriz con \((i,j)\) elemento \(\gamma_{i,j}(t)\)

Una matriz finita de Markov satisface la ecuacion Chapman-Kolmogorov: \[\Gamma(t+u)=\Gamma(t)+\Gamma(u)\] La ecuancion de Chapman-kolmogorov implica para todo \(t\in\mathbb{N}\), \[\Gamma(t)=\Gamma(1)^t\] sea esto la matriz de probabilidad de transicion de t-pasos es la t-esima potencia de \(\Gamma(1)\), la matriz de primer paso. La matriz \(\Gamma(1)\), es la abreviacion de \(\Gamma\) es una matriz cuadrada cuyas filas suman 1. \[\begin{pmatrix} \gamma_{11} & \dots & \gamma_{1n} \\ \vdots & \ddots & \vdots \\ \gamma_{m1} & \dots & \gamma_{mn}\\ \end{pmatrix}\] Donde \(m\) denota el numero de estados de la cadena de Markov. La sentencia de las filas sumando 1 puede ser escrito como \(\Gamma1' = 1'\) esto es, el vector columna \(1'\) es un eigenvector de \(\Gamma\) y corresponden a los eigenvalores 1. Nos referimos a \(\Gamma\) como la matriz de probabilides de transicion (t.p.m, por sus siglas en ingles)

La probabilidad no condicionada \(Pc(C_t=j)\) de una cadena de Markov, en un estado dado, en un momento t es de interes y es denotado por el vector fila:

\(u(t)=(Pr(C_t=1),...,Pr(C_t=m)), \in \mathbb{N}\) Es decir nos permite saber el estado actual de la matriz. Siendo \(u(1)\) la distribucion incial de la matriz, para deducir la distribucion en el tiempo \(t+1\) desde \(t\) solo multiplicamos por la matriz de probabilidad de transicion: \[u(t+1)=u(1)\Gamma\]

Distribuciones Estacionarias

Una cadena de Markov con matriz de probabilidad de transicion \(\Gamma\) se dice que tiene una distribucion \(\delta\) (un vector fila con elementos no negativos) si \(\delta\Gamma=\delta\) y \(\delta1'=1\).

El primer requerimiento expresa la estacionalidad y el segundo que \(\delta\) es de echo una distribucion de probabilidad.

La matriz \[\begin{pmatrix} 1/3 & 1/3 & 1/3 \\ 2/2 & 0 & 1/3 \\ 1/2 & 1/2 & 0 \\ \end{pmatrix}\]

Tiene una distribucion estacionaria \(\delta=\frac{1}{32}(15,9,8)\). Siendo que \(u(t+1)=u(t)\Gamma\) una cadena de Markov empezando desde esta distribucion tendra la misma distribucion en cualquier punto en el tiempo, esto es una matriz de Markov estacionaria.

Una cadena de Markov irreducible (homogenea, discreta en tiempo, espacio de estados finitos) tiene una unica, estrictamente positiva distribucion estacionaria. El vector \(\delta\) con elementos no negativos es una distribucion estacionaria de la cadena de Markov con t.p.m \(\Gamma\), si y solo si: \[\gamma(I_m-\Gamma +U)=1\] Donde 1 es un vector de unos, \(I_m\) es la matriz identidad \(m\times m\) y \(U\) es la matriz \(m\times m\) de unos.

La matriz estacionaria puede ser encontrada removiendo una de las ecuaciones del sistema \(\delta\Gamma=\delta\) y reemplazandola por \(\sum_i\delta_i=1\)

Funcion de Autorrelacion

Asumimos que los estados \(1,2,..,m\) son cuantitativos y no categoricos, la funcion de autocorrelacion de \({C_t}\), se asume que es estacionaria e irreducible, se puede obtener que:

Primero, v=(1,2,…,\(m\)) y V=diag(1,2,…,\(m\)) tenemos los enteros no negativos \(k\), \[Cov(C_t,C_t+k)0\delta V\Gamma^kv'-(\delta v')^2\]

Si \(\Gamma\) es diagonalizable, y los eigenvalores (diferentes a 1) se denotan por \(\omega_2,\omega_3,...,\omega_m\), entonces \(\Gamma\) puede ser escrito como: \[\Gamma=U\Omega U^-1\] Donde \(\Omega\) es la diag(\(1,\omega_2,\omega_23,...,\omega_m\)) y U corresponde a los eigenvectores de \(\Gamma\). Tenemos para los valores no negativos \(k\) \[Cov(C_t,C_{t+k})=\delta VU\Omega ^kU^{-1}v'-(\delta v')^2\] \[= a\Omega ^kb'-a_1b_1\] \[\sum_{i=2}^{m}a_ib_i\omega_i^k\] donde \(a=\delta VU\) y \(b'=U^{-1}v'\) \[\rho(k)\equiv Corr(C_t,C_{t+k})=\frac{\sum_{i=2}^{m}a_ib_i\omega_i^k}{\sum_{i=2}^{m}a_ib_i}\]

Cadenas de Markov

Una vez observamos una cadenas de Markov y queremos estimar las probabilidades de transicion, una solucion - pero no la unica - is encontrar las transiciones y estimar las probabilidades como frecuencias relativas.

Suponiendo entonces, que queremos estimar \(m^2-m\) parametros \(\gamma_{ij} (i\neq j)\) de un estado \(m\) de la cadena de Markov \({C_t}\) El “likelihood” condicionado de la primera observacion es: \[L=\prod_{i=1}^{m}\prod_{j=1}^{m}\gamma_{ij}^{f_{ij}}\] Siendo el Log de este: \[l=\sum_{i=1}^{m}\sum_{j=1}^{m}f_{ij} log \gamma_{ij}=\sum_{i=1}^{m}l_i\]

y maximizamos \(l\), maximizando cada \(l_i\) por separado. Sustituyendo \(1-\sum_{k\neq i} \gamma_{ik}\), diferenciando \(l_i\) respecto a la diagonal de la probabilidad de transición \(\gamma_{ij}\) e igualando a cero tenemos que: \[0=\frac{-f_{if}}{1-\sum_{k\neq i}\gamma_{ik}}+\frac{f_{ij}}{\gamma_{ij}}=\frac{f_{ii}}{\gamma_ii}+\frac{f_{ij}}{\gamma_{ij}}\] A menos que un denominador sea cero, \(f_{ij}\gamma_{ii}=f_{ii}\gamma_{ij}\) y \(\sum_{j=1}^{m}f_{ij}=f_{ii}\). Esto implica que en el maximo del “likelihood”, \[\gamma_{ii}=\frac{f_{ii}}{\sum_{j=1}^{m}f_{ij}}\text{ y } \gamma_{ij}=\frac{f_{ij}\gamma_{ii}}{f_{ii}}=\frac{f_{ii}}{\sum_{i=1}^{}m}f_ij\] El estimador \(\gamma=\frac{f_{in}}{\sum_{k=1}^{m}f_{ik}(i,j=1,...,m)}\), es el la probabilidad de transicion empirica, el estimador \(\Gamma\) satisface el requerimiento de la sumas de las filas se igual a 1

Suponiendo la cadena de Markov \({C_t}\) toma valores 0 y 1, y deseasmos estimar las probabilidades de transicion \(\gamma_{ij}\) de la secuencia de observaciones en que \(f_{ij}\) transiciones del estado \(i\) al estado \(j\) (\(i\),\(j\)=0,1) y \(f_{11}>0\) pero \(f_{00}=0\). entonces en las observaciones un cero es seguido por un uno. \(c=f_{10}+(1-C_1)\) y \(d=f_{11}\). Las probabilidades de transicions estan dadas por \[\gamma_{01}=1 \text{ y } \gamma_{10}=\frac{-(1+d)+((1+d)^2 +4c(c+d-1))^{\frac{1}{4}}}{2(c+d-1)}\]

Cadenas de Markov de orden superior

En el caso donde las observaciones con estados finitos parecen no satisfacer las propiedades de Markov, una alternativa es utiliza cadenas de Markov de orden superior, satisfaciendo: \[Pr(C_t|C_{t-1},C_{t-2,...})=Pr(C_t|C_{t-1},...,C_{t-1})\] Aunque el modelo no es una cadena de Markov per se, podemos redefinir el modelo de forma que el resultante lo sea. Si dejamos \(Y_t=(C_{t-l+1},C_{t-l+2},...,C_t)\) entonces \({Y_t}\) es una cadena de Markov de primer orden en \(M^t\) en el espacio de estados de \({C_t}\).

Una cadena de Markov de segundo orden, es estacionaria si se caracteriza por las probabilidades de transicion: \[\gamma(i,l,k)=Pr(C-t=k|C_{t-1},C_{t_2}=i)\] y tiene distribucion estacionaria bivariable \(u(j,k)=Pr(C_{t-1}=j,C_t=k)\) satisfaciendo \[u(j,k)=\sum_{i=1}^{m}u(i,j)\gamma(i,j,k) \text{ y } \sum_{j=1}^{m}\sum_{k=1}^{m}u(j,k)=1\] El uso general de las cadenas de Markov de orden superior incrementa el numero e parametros del modelo; una cadena general de Markov de orden \(l\) en \(m\) estados tiene \(m^l(m-1)\) probabilidades de transicion independientes. Param \(m=2\) hay modelos equivalentes, pero \(M>2\) pueden ser representados en una variedad de patrones dependientes y estructuras de autocorrelacion. En ambos casos un incremento de 1 orden en la cadena de Markov requiere solo un parametro adicional.

Los modelos de Raftery, quien denomina modelos “distribuciones de transicion mixtas” estan definidos como los procesos \({C_t}\) que toman valores \(M={1,2,...,m}\) y satisfacen \[Pr(C_t=j_0|C_{t-1}=j_1,...,C_{t-l}=jl) = \sum_{i=1}^{l}\lambda_i q(j_i,j_0)\] Donde \(\sum_{i=1}{l}\lambda=1\) y \(Q=(q(j,k))\) es una matriz \(m\times m\) con elementos no negativos que sus filas suman 1; cuyo lado derecho de la ecuacion esta acotada entre ceo y uno para todos los \(j_0,j_1,...,j_l\in M\)

CAPITULO 2

Modelo de Markov escondidas

Definiciones y propiedades

Modelo simple de Markov escondido

Considerando los visto en la serie de terremotos en el capitulo anterior, se puede apreciar que las observaciones no tienen un limite de cantidad. Debido a esto, utilizar la distribucion de Poisson es la decision natural. En el capitulo uno se considero como puede afectar la interpretacion de la serie como una mezcla de distribuciones de Poisson, Cada una con una media \(\lambda_m\) donde m se refiere a la serie de datos y cada distribucion tiene un posibilidad de \(\delta_m\) donde \[\sum_{ i = 1}^{\infty} \delta_i = 1\]

Un modelo como una mezcla independiente se considera que no representara correctamente el modelo de los terremotos debido a que por definicion estas mezcla no tienen dependecia alguna con los valores anteriores sin embargo de acuerdo con los resultados de la grafica ACF una dependencia muy baja entre las observaciones, pero por fines ilutrativos se obtendra el valor de la correlaciones.

muestra <- read_xlsx("EarthQuakeData.xlsx", sheet="Sheet2", skip = 2)[-109,] %>% 
  rename("x" = !!names(.[2])) %>% 
  select(x)

acf(muestra$x)

Lo Basico

La cadena escondida de Markov es una mezcla particular debido a su dependencia. Donde toda observacion \(X_t\) tiene adicionalmente una categoria \(C_t\) la cual se indica a cual de las distribuciones pertenece mientras que se puede calcular las probabilidades de la siguiente manera: \[Pr(C_t|C_{t-1}), t = 2,3..\] \[Pr(X_t|X_{t-1}, C_{}) t = 2,3.\]

Para este modelo se deben considerar dos parametros, primero un parametro no observado del proceso, \(C_t\) que satisface todas las propiedades de Markov, y un segundo parametro dependiente del estado, \(X_t\) en el cual la distribución de \(X\) depende unicamente de estado actual, \(C\).

Si la cadena de Markov tiene m estados, entonces se le conoce como HMM de m-estados. En este caso estamos considerando que apesar que ya se calcularos los valores de \(\delta_i\), existe adicionalmente una matriz de transicion como: \[ \Gamma = \begin{bmatrix} \gamma_{11} & \gamma_{12} \\ \gamma_{21} & \gamma_{22} \\ \end{bmatrix}\]

En contraste, en un caso de mezclas independientes, la distribucion de \(C_t\), el estado en el tiempo \(t\), si depende de \(C_{t-1}\). Por lo tanto, para considerar como se pueden representar las distribuciones de maneras discretas y continuas se debe puede utilizar como base la siguiente expresion para el estado \(i = 1,2,3..\) \[p_i(x)= Pr(X_t=x | C_t = i) \]

\(p_i\) es la probabilidad de \(X_t\) si la cadena de Markov esta en el estado \(i\) en el momento \(t\). En caso de variables continuas es necesarios considerar que la probabilidad como una funcion de densidad. A estas distribuciones \(p_i\) se les consideran como distribuciones dependientes del estado, teniendose \(m\) estados.

Distribucion Marginal

Ocacionalmente se necesita la distribucion marginal de \(X_t\) y distribuciones marginales mas elevadas, tal como \((X_t, X_{t+k})\). Estos resultados deben ser derivados de los casos de la matrices homogenias pero que no son necesariamente estacionarias.

Univariate distributions

Para valores discretos observados \(X_t\), definiendo \(u_i(t) = Pr(C_t = i)\) para \(t = 1,2 ...T\) se tiene

\[Pr(X_t = x) = \sum_{i = 1}^{m} Pr(C_t = i) Pr(X_t = x| C_t = i) = \sum_{i = 1}^{m} u_i(t)p_i(x) \]

La expreción se puede escribir de forma matricial de la siguiente manera: [Pr(X_t = x) = (u_1(t)…u_m(t)) \[\begin{bmatrix} p_1(x) & & 0 \\ & \ddots & \\ 0 & & p_m(x) \\ \end{bmatrix}\] \[\begin{bmatrix} 1 \\ \vdots \\ 1 \\ \end{bmatrix}\]

]

Donde \(P(x)\) es definida como la diagonal de los \(p_i(x)\). Considerando que \(u_i(t) = u(1)\Gamma^{t-1}\), se obtiene: \[Pr(X_t = x) = u(1)\Gamma^{t-1}P(x) 1^{`}\]

Pero si se consideran el cambio de estado como una distribución estacionaria entonces se puede reescribir teniendo en cuenta que: \(\delta\Gamma^{t-1} = \delta\) para todo \(t \in \mathbb{N}\) se obtiene:

\[Pr(X_t = x) = \delta P(x)1^{`}\]

Momentos

Primero debemos observar que: \[E(X_t) = \sum_{i=1}^{m}E(X_t|C_t=i)Pr(C_t=i)\]

Que al considerarla como un caso estacionario seria: \[E(X_t) = \sum_{i=1}^{m}\delta_i E(X_t|C_t=i)\]

En base a la ecuacion de arriba y considerado las propiedades de los momentos y de las cadenas se obtienen los siguientes momentos:

\[E(X_t) = \delta_1 \lambda_1 + \delta_2 \lambda_2\] \[Var(X_t) = E(X_t) \delta_1 \delta_2 ( \lambda_1 +\lambda_2)^2 > E(X_t)\] \[Cov(X_t, X_{t+k}) = \delta_1 \delta_2 ( \lambda_1 +\lambda_2)^2(1-\gamma_{12}-\gamma_{21})^k , For k \in \mathbb{N} \]

Cabe mencionar que la formula resultante de la covarianza tiene una forma tal que cuando las medias son iguales \(\lambda_1 = \lambda_2\) la covarianza se vuelve cero.

Likelihood

El objetivo de la seccion desarrollar la formula del likelihood \(L_T\) de \(T\) observaciones consecuentes asumiendo que son el resultado de una HMM de m estados, Esta formula existe pero para obtenerse el resultado es necesario calcular al rededor de \(m^T\) terminimos de \(2T\) factores, aparentando requerir \(O(Tm^T)\). Sin embargo hace tiempo esta comprobado que se puede calcular con \(O(Tm^2)\) operaciones. Esto es posible mediante la estimacion de los parametros por medio de optimizacion numerica.

Likelihood en general

Consdierando que el likelihood de una HMM en general. Supongamos una secuencias de observaciones \(x_1, x_2, .... x_T\) generador por este tipo de modelo. Se busca la probabilidad \(L_T\) de observar esta secuencia, como calculado bajo el modelo m-estado HMM que inicialmente tiene una distribución \(\delta\) y \(\Gamma\) de la cadena de Markov.

CAPITULO 3

Estimacion por maximizacion direta del likelihood

El likelihood de una Cadena de Markov Escondida (HMM), se puede representar por la siguiente ecuacion: \[L_t = P_r (X^{(T)} = x^{(T)}) = \delta P(x_1) \Gamma P(x_2) ... \Gamma P(x_T)1'\]

Siendo \(\delta\) la distribución inicial y \(P(x)\) la matriz diagonal de \(m*m\) con densidad de probabilidad \(p_i(x)\)

\(L_T\) se puede calcular como \(L_T = \alpha_T1'\) tomando \(\alpha_1 = \delta P(x_1)\) y de forma recursiva:

\(\alpha_t = \alpha_{t-1} \Gamma P(x_t),\) para \(t=2,3,...,T\)

Si la cadena de Markov se asume como estacionaria \((\delta=\delta \Gamma)\), en ese caso \(\alpha_0 = \delta\), por lo tanto:

\(\alpha_t = \alpha_{t-1} \Gamma P(x_t),\) para \(t=1,2,...,T\)

Escalacion del calculo del likelihood

El likelihood es un producto de matrices, no de escalares, es por ello que no se puede simplemente calcular el log del likelihood como la suma de los logs de todos sus factores. Se sugiere un metodo de calcular el likelihood que supone \(log(p+q)\), en donde \(p>q\),

\[log(p+q) = log(p)+log(1 + q/p) = log(p)+log(1+exp(log(q)-log(p)))\]

El logaritmo de \(L_T\) se calcula escalando el vector de probabilidades \(\alpha_t\), es por ello que el vector se escala en cada \(t\) entonces sus elementos se suman a 1, manteniendo la sumatoria de los logs de los factores escalados. Por lo tanto el \(log(L_T)\) se puede definir como:

\[log(L_T) = \sum_{t=1}^T log(w_t/w_{t-1}) = \sum_{t=1}^T log(\phi_{t-1}\Gamma P(x_t)1')\]

Para calcular el log-likelihood de las observaciones \(x_1, x_2,...,x_t\) dentro del modelo Poisson-HMM con, al menos 2 estados, con matriz de probabilidades de transicion \(\Gamma\), vector de medias de estados dependientes \(\lambda\), y distribucion inicial \(\delta\), se utiliza el siguiente codigo implementado en R.

alpha <- params$delta * dpois(x[1], params$lambda)
lscale <- log(sum(alpha))
alpha <- alpha/sum(alpha)

for (i in 2:T) {
  #alpha <- alpha %*% Gamma*as.numeric(dpois(x[i], params$lambda))
  lscale <- lscale+log(sum(alpha))
  alpha <- alpha/sum(alpha)
}
lscale
## [1] -2.386378

Maximizacion del likelihood sujeto a restricciones

Primero se debe realizar una reparametrizacion para evitar las restricciones. Los elementos de \(\Gamma\), los de \(\lambda\) y el vector de medias de estados dependientes en una Poisson-HMM estan sujetos a restricciones de no negatividad y de otro tipo.

Las restricciones se pueden clasificar en 2 grupos. El primer grupo de restricciones son aquellas que dependen del estado de la distrubucion que se haya elegido, por ejemplo, la probabilidad de exito de una distribucion binomial se encuentra entre 0 y 1.

En el caso Poisson-HMM las restricciones mas importantes son: 1. La media \(\lambda_i\) del estado dependiente de la distribucion debe ser no negativa, para \(i=1,...,m\) 2. Las filas de la matriz de probabilidades de transicion \(\Gamma\) debe agregarse a 1, y todos los parametros \(\gamma_{ij}\) deben ser no negativos.

Para hacer la transformacion de los parametros, vamos a definir \(\eta = log \lambda_i\), para \(i=1,...,m\) en donde \(\eta \in \mathbb{R}\)

Luego de obtener la maximizacion del likelihood sin restricciones, se pueden obtener los parametros con restricciones transformando nuevamente: \(\lambda_i = exp (\eta_i)\).

La reparametrizacion de la matriz \(\Gamma\) se hace de la siguiente manera: \[\sum_{j=1}^m \gamma_{ij} = 1 \hspace{1cm} (i=1,...,m)\]

\(\Gamma\) tiene \(m^2\) entradas pero unicamente \(m(m-1)\) parametros libres, por lo que las se puede transformar a numeros reales sin restricciones \(\tau_{ij}, i \neq j\).

Para el caso \(m=3\) se define la matriz \[T = \begin{pmatrix} - & \tau_{12} & \tau_{13} \\ \tau_{21} & - & \tau_{23} \\ \tau_{31} & \tau_{32} & - \\ \end{pmatrix}\]

en donde \(\tau{ij} \in \mathbb{R}\). De tal forma que \(g : \mathbb{R} \rightarrow \mathbb{R}^+\) sea estrictamente una funcion creciente, \[g(x) = e^x \hspace{0.5cm} o \hspace{0.5cm} g(x) = \Bigg \{ \begin{array}{} e^x & x \leq0 \\x+1 & x \geq0 \end{array}\]

Entonces: \[ v_{ij} = \Bigg \{ \begin{array}{} g(\tau_{ij}) & para & i \neq j \\1 & para & i = j \end{array}\]

y:

\(\gamma_{ij} = \frac {v_{ij}} {\sum_{k=1}^m v_{ik}} \hspace{0.5cm} para \hspace{0.5cm} i,j = 1,2,...,m\) y \(\Gamma = (\gamma_{ij})\)

Por lo tanto los parametros \(\eta_i\) y \(\tau_{ij}\) los llamaremos parametros de trabajo, mientras que los paramentros \(\lambda_i\) y \(\gamma_{ij}\) seran los parametros naturales.

Utilizando estas transformaciones para \(\Gamma\) y \(\lambda\) el calculo de los parametros de la maximizacion del likelihood se pueden calcular asi: 1. Maximizar \(L_T\) con respecto a los parametros de trabajo \(T = \{\tau_{ij}\}\) y \(\eta = (\eta_1,...,\eta_m)\) 2. Transformar los estimadores de los parametros de trabajo en los estimadores de los parametros naturales: \(T \rightarrow \Gamma\), \(\eta \rightarrow \lambda\).

A continuacion se presenta un ejemplo de la una funcion en codigo de R haciendo referencia a la transformacion de los parametros naturales de Poisson en parametros de trabajo y viceversa.

Transformacion de los parametros naturales de Poisson en parametros de trabajo

pois.HMM.pn2pw <- function(m,lambda,gamma,delta=NULL,stationary=TRUE)
{
 tlambda <- log(lambda)
 foo     <- log(gamma/diag(gamma))
 tgamma  <- as.vector(foo[!diag(m)])
 if(stationary) {tdelta <-NULL} else {tdelta<-log(delta[-1]/delta[1])}
 parvect <- c(tlambda,tgamma,tdelta)
 return(parvect)
}

Transformacion de los parametros de trabajo de Poisson en parametros naturales

pois.HMM.pw2pn <- function(m,parvect,stationary=TRUE)
{
 lambda        <- exp(parvect[1:m])
 gamma         <- diag(m)
 gamma[!gamma] <- exp(parvect[(m+1):(m*m)])
 gamma         <- gamma/apply(gamma,1,sum)
 if(stationary) {delta<-solve(t(diag(m)-gamma+1),rep(1,m))} else
                {foo<-c(1,exp(parvect[(m*m+1):(m*m+m-1)]))
                delta<-foo/sum(foo)}
 return(list(lambda=lambda,gamma=gamma,delta=delta))
}